home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / SPCTRM.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  69 lines
  1. PROCEDURE spctrm(VAR p: glmarray; m,k: integer; ovrlap: boolean;
  2.        VAR w1: gl4marray; VAR w2: glmarray);
  3. (* Programs using routine SPCTRM must include 'dfile' as
  4. a main routine parameter, and declare
  5. VAR
  6.    dfile: text;
  7. They must open and close this file appropriately.
  8. Also they must define data types
  9. TYPE
  10.    glmarray = ARRAY [1..m] OF real;
  11.    gl4marray = ARRAY [1..4*m] OF real;
  12. and must declare two workspace matrices w1,w2 with types as shown
  13. in the arguments above. *)
  14. VAR
  15.    mm,m44,m43,m4,kk,joffn,joff,j2,j,jj: integer;
  16.    w,sumw,facp,facm,den: real;
  17. FUNCTION window(j: integer; facm,facp: real): real;
  18.    BEGIN
  19.       window := (1.0-abs(((j-1)-facm)*facp))      (* Parzen *)
  20. (*      window := 1.0            *)   (* Square *)
  21. (*      window := (1.0-sqr(((j-1)-facm)*facp))    *)   (* Welch  *)
  22.    END;
  23. BEGIN
  24.    mm := m+m;
  25.    m4 := mm+mm;
  26.    m44 := m4+4;
  27.    m43 := m4+3;
  28.    den := 0.0;
  29.    facm := m-0.5;
  30.    facp := 1.0/(m+0.5);
  31.    sumw := 0.0;
  32.    FOR j := 1 TO mm DO sumw := sumw+sqr(window(j,facm,facp));
  33.    FOR j := 1 TO m DO p[j] := 0.0;
  34.    IF (ovrlap) THEN BEGIN
  35.       FOR j := 1 TO m DO read(dfile,w2[j])
  36.    END;
  37.    FOR kk := 1 TO k DO BEGIN
  38.       FOR joff := -1 TO 0 DO BEGIN
  39.          IF (ovrlap) THEN BEGIN
  40.             FOR j := 1 TO m DO w1[joff+j+j] := w2[j];
  41.             FOR j := 1 TO m DO read(dfile,w2[j]);
  42.             joffn := joff+mm;
  43.             FOR j := 1 TO m DO w1[joffn+j+j] := w2[j] END
  44.          ELSE BEGIN
  45.             FOR jj := 0 TO ((m4-joff-2) DIV 2) DO BEGIN
  46.                j := joff+2+2*jj;
  47.                read(dfile,w1[j])
  48.             END
  49.          END
  50.       END;
  51.       FOR j := 1 TO mm DO BEGIN
  52.          j2 := j+j;
  53.          w := window(j,facm,facp);
  54.          w1[j2] := w1[j2]*w;
  55.          w1[j2-1] := w1[j2-1]*w
  56.       END;
  57.       four1(w1,mm,1);
  58.       p[1] := p[1]+sqr(w1[1])+sqr(w1[2]);
  59.       FOR j := 2 TO m DO BEGIN
  60.          j2 := j+j;
  61.          p[j] := p[j]+sqr(w1[j2])+sqr(w1[j2-1])
  62.             +sqr(w1[m44-j2])+sqr(w1[m43-j2])
  63.       END;
  64.       den := den+sumw
  65.    END;
  66.    den := m4*den;
  67.    FOR j := 1 TO m DO p[j] := p[j]/den
  68. END;
  69.